Location-scale Meta-analysis and Meta-regression as a Tool to Capture Large-scale Changes in Biological and MethodologicalHeterogeneity: a Spotlight on Heteroscedasticity

1 Introduction (Outline)

Here, we illustrates how to fit location-scale meta-analysis and meta-regression models in R, using brms (Bayesian approach), demonstrating how to model both mean (location) and variance (scale) components under a meta-analytic framework.

We cover: 1. A dataset with a categorical moderator (e.g., habitat type). 2. A dataset with a continuous moderator (e.g., log-elevation). 3. An example showing how to examine publication biases in both location and scale parts (e.g., small-study divergence, Proteus effect).

For each example, we start with fitting the following location-scale model (or double-hierarchical model):

\[ y_i = \beta_0^{(l)} + u_{j[i]}^{(l)} + e_i^{(l)} + m_i^{(l)}, \] \[ \ln(\sigma_{e_i}) = \beta_0^{(s)} + u_{j[i]}^{(s)}. \] \[ \begin{equation} \begin{pmatrix} u_j^{(l)} \\ u_j^{(s)} \end{pmatrix} \sim \mathcal{N}\left( \begin{pmatrix}0 \\ 0\end{pmatrix}, \begin{pmatrix} \sigma_{u(l)}^{2} & \rho_u \sigma_{u(l)} \sigma_{u(s)} \\ \rho_u \sigma_{u(l)} \sigma_{u(s)} & \sigma_{u(s)}^{2} \end{pmatrix} \right). \end{equation} \]

\[ e_i \sim \mathcal{N}(0,\sigma_e^2), \space \text{and} \space m_i \sim \mathcal{N}(0,\sigma_{m_i}^2). \]

TODO - description….

For Example 3, we also use:

\[ \begin{equation} \begin{pmatrix} u_j^{(l)} \\ u_j^{(s)} \end{pmatrix} \sim \mathcal{N}\left( \begin{pmatrix}0 \\ 0\end{pmatrix}, \begin{pmatrix} \sigma_{u(l)}^{2} & 0 \\ 0 & \sigma_{u(s)}^{2} \end{pmatrix} \right). \end{equation} \]

Then, we fit the following location-scale meta-regression

\[ y_i = \beta_0^{(l)} + \beta_1^{(l)} x_{1i} + \cdots + \beta_p^{(l)} x_{pi} + u_{j[i]}^{(l)} + e_i^{(l)} + m_i^{(l)}, \] \[ \ln(\sigma_{e_i}) = \beta_0^{(s)} + \beta_1^{(s)} x_{1i} + \cdots + \beta_p^{(s)} x_{pi}, \]

\[ u_j^{(l)} \sim \mathcal{N}(0,\sigma_{u(l)}^{2}), \space e_i \sim \mathcal{N}(0,\sigma_e^2), \space \text{and} \space m_i \sim \mathcal{N}(0,\sigma_{m_i}^2) \] TODO - description….

For Example 1, we also fit a phylogenetic location-scale meta-regression as follows:

\[ y_i = \beta_0 + a_{k[i]} + s_{k[i]} + u_i + e_i + m_i \]

\[ \ln(\sigma_{e_i}) = \beta_0^{(s)} + \beta_1^{(s)} x_{1i} + \cdots + \beta_p^{(s)} x_{pi}, \] \[ a_k^{(l)} \sim \mathcal{N}(0,\sigma_{a(l)}^{2}\mathbf{A}), \space \text{and} \space s_k^{(l)} \sim \mathcal{N}(0,\sigma_{s(l)}^{2}). \]

with

\[ u_j^{(l)} \sim \mathcal{N}(0,\sigma_{u(l)}^{2}), \space e_i \sim \mathcal{N}(0,\sigma_e^2), \space \text{and} \space m_i \sim \mathcal{N}(0,\sigma_{m_i}^2). \]
TODO - description….

Later, we use metafor (a frequentist implementation) and blsmeta (another Bayesian implementation) to check if results from brms match those from metafor and blsmeta when we use the same model.

Please refer to the associated paper for theoretical background, formulas, and details on location-scale models in ecology and evolution.

2 Prerequisites

2.1 Loading pacakges

# Attempt to load or install necessary packages
if(!require(pacman)) install.packages("pacman")
pacman::p_load(
  tidyverse,
  tidybayes,
  here,
  patchwork,
  orchaRd,  # see the note
  ggplot2,
  pander,
  brms,
  metafor,
  blsmeta # see the note
)

Note

not about functions to install…..

It’s important…

2.2 Adding custum functions

These functions are to visualize brms results

# Function to get variable names dynamically
get_variables_dynamic <- function(model, pattern) {
  variables <- get_variables(model)
  variables[grep(pattern, variables)]
}

rename_vars <- function(variable) {
  variable <- gsub("b_Intercept", "b_l_int", variable)
  variable <- gsub("b_sigma_Intercept", "b_s_int", variable)  
  variable <- gsub("b_habitatterrestrial", "b_l_contrast", variable)            
  variable <- gsub("b_sigma_habitatterrestrial", "b_s_contrast", variable)
  variable <- gsub("b_methodpersisitent", "b_l_contrast", variable)            
  variable <- gsub("b_sigma_methodpersisitent", "b_s_contrast", variable)
  variable <- gsub("b_elevation_log", "b_l_slope", variable)            
  variable <- gsub("b_sigma_elevation_log", "b_s_slope", variable)
  variable <- gsub("sd_study_ID__Intercept", "sd_l_study_ID", variable)            
  variable <- gsub("sd_study_ID__sigma_Intercept", "sd_s_study_ID", variable)
  variable <- gsub("cor_study_ID__Intercept__sigma_Intercept", "cor_study_ID", variable)  
  return(variable)
}

# Function to visualize fixed effects
visualize_fixed_effects <- function(model) {
  fixed_effect_vars <- get_variables_dynamic(model, "^b_")
  if (length(fixed_effect_vars) == 0) {
    message("No fixed effects found")
    return(NULL)
  }
  
  tryCatch({
    fixed_effects_samples <- model %>%
      spread_draws(!!!syms(fixed_effect_vars)) %>%
      pivot_longer(cols = all_of(fixed_effect_vars), names_to = ".variable", values_to = ".value") %>%
      mutate(.variable = rename_vars(.variable))
    
    ggplot(fixed_effects_samples, aes(x = .value, y = .variable)) +
      stat_halfeye(
        normalize = "xy", 
        point_interval = "mean_qi", 
        fill = "lightcyan3", 
        color = "lightcyan4"
      ) +
      geom_vline(xintercept = 0, linetype = "dashed", color = "#005") +
      labs(y = "Fixed effects", x = "Posterior values") +
      theme_classic()
  }, error = function(e) {
    message("Error in visualize_fixed_effects: ", e$message)
    return(NULL)
  })
}

# Function to visualize random effects
visualize_random_effects <- function(model) {
  random_effect_vars <- get_variables_dynamic(model, "^sd_")
  random_effect_vars <- random_effect_vars[random_effect_vars != "sd_es_ID__Intercept"]
  if (length(random_effect_vars) == 0) {
    message("No random effects found")
    return(NULL)
  }
  
  tryCatch({
    random_effects_samples <- model %>%
      spread_draws(!!!syms(random_effect_vars)) %>%
      pivot_longer(cols = all_of(random_effect_vars), names_to = ".variable", values_to = ".value") %>%
      mutate(.variable = rename_vars(.variable)) #%>%
      #mutate(.value = .value)  # leave SD as it is
    
    ggplot(random_effects_samples, aes(x = .value, y = .variable)) +
      stat_halfeye(
        normalize = "xy", 
        point_interval = "mean_qi", 
        fill = "olivedrab3", 
        color = "olivedrab4"
      ) +
      geom_vline(xintercept = 0, linetype = "dashed", color = "#005") +
      labs(y = "Random effects (SD)", x = "Posterior values") +
      theme_classic()
  }, error = function(e) {
    message("Error in visualize_random_effects: ", e$message)
    return(NULL)
  })
}

# Function to visualize correlations
visualize_correlations <- function(model) {
  correlation_vars <- get_variables_dynamic(model, "^cor_")
  if (length(correlation_vars) == 0) {
    message("No correlations found")
    return(NULL)
  }
  
  tryCatch({
    correlation_samples <- model %>%
      spread_draws(!!!syms(correlation_vars)) %>%
      pivot_longer(cols = all_of(correlation_vars), names_to = ".variable", values_to = ".value") %>%
      mutate(.variable = rename_vars(.variable))
    
    ggplot(correlation_samples, aes(x = .value, y = .variable)) +
      stat_halfeye(
        normalize = "xy", 
        fill = "#FF6347", 
        color = "#8B3626"
      ) +
      geom_vline(xintercept = 0, linetype = "dashed", color = "#005") +
      labs(y = "Correlations", x = "Posterior values") +
      theme_classic()
  }, error = function(e) {
    message("Error in visualize_correlations: ", e$message)
    return(NULL)
  })
}

3 Example 1: Categorical Moderator (Thermal Tolerance Dataset)

Data source (illustrative): Pottier, P. et al. (2022). Developmental plasticity in thermal tolerance: Ontogenetic variation, persistence, and future directions. Ecology Letters, 25(10), 2245–2268.

3.1 Loading the data

This dataset (thermal.csv) has an effect size dARR (the developmental acclimation response ratio: the ratio of the slopes of acclimation) and sampling variance Var_dARR. Key moderators are two categorical variables: (1) habitat (e.g., aquatic vs. terrestrial: where animals live) and (2) method (initial vs. persisitent: experimental design where only temparature increase was applied during intial phaes or over the entire experimental perid).

dat <- read.csv(here("data", "thermal.csv"))

# selecting varaibles we need
dat <- dat %>% select(dARR, Var_dARR, es_ID, population_ID, study_ID, exp_design, habitat)

# creating SE (sqrt(Var))
dat$si <- sqrt(dat$Var_dARR)

# creating a new variable (method) according to the oringal paper
dat$method <- ifelse(dat$exp_design == c("A", "B", "C"), "initial", "persisitent")

#head(dat)
str(dat)
## 'data.frame':    1089 obs. of  9 variables:
##  $ dARR         : num  0.0528 0.0428 0.054 0.0115 0.0503 ...
##  $ Var_dARR     : num  0.000171 0.000192 0.000239 0.000131 0.000175 ...
##  $ es_ID        : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ population_ID: int  1 1 2 2 2 3 3 3 3 3 ...
##  $ study_ID     : int  1 1 1 1 1 2 2 2 2 2 ...
##  $ exp_design   : chr  "D" "D" "D" "D" ...
##  $ habitat      : chr  "terrestrial" "terrestrial" "terrestrial" "terrestrial" ...
##  $ si           : num  0.0131 0.0139 0.0154 0.0114 0.0132 ...
##  $ method       : chr  "persisitent" "persisitent" "persisitent" "persisitent" ...

# If needed, create effect-size IDs
# dat$es_ID <- as.factor(1:nrow(dat))
# dat$study_ID <- <some factor> if multiple effect sizes per study

3.2 Creating a variance-covariance matrix for sampling errors

3.2.1 Assuming indepedence among effect sizes

If each effect size has an independent sampling error, we can build a diagonal matrix:

vcv <- diag(dat$Var_dARR)
rownames(vcv) <- dat$es_ID
colnames(vcv) <- dat$es_ID

3.2.2 Modeling non-independence

In reality, non-independence of sampling errors are common (effect sizes are obtained partially or all from the same subjects or individuals)

# here correlated structure or cluster is population_ID
# vcalc is from the metafor pacakge
# we will not sue this for models but this is for demonstration

vcv2 <- vcalc(vi = Var_dARR, 
             cluster = population_ID, 
             obs = es_ID, rho = 0.5, 
             data = dat)
rownames(vcv2) <- dat$es_ID
colnames(vcv2) <- dat$es_ID

3.3 Multilevel location-scale meta-analysis

#TODO - may need to rerun this model

# p in the random effects in location and scale parts allows correlation to be modeled 
form_ma1 <- bf(dARR
            ~ 1   +
              (1|p|study_ID) + # this is u_l (the between-study effect)
              (1|gr(es_ID, cov = vcv)), # this is m (sampling error)
            sigma ~ 1 + 
              (1|p|study_ID) # residual = the within-study effect goes to the scale
)


# Generate default priors
prior_ma1 <- default_prior(form_ma1, 
                          data=dat, 
                          data2=list(vcv=vcv),
                          family=gaussian())
prior_ma1$prior[5] = "constant(1)" # meta-analysis assumes sampling variance is known so fixing this to 1
prior_ma1

# fitting model
fit_ma1 <- brm(
  formula = form_ma1,
  data = dat,
  data2 = list(vcv=vcv),
  chains = 2,
  cores = 2,
  iter = 6000,
  warmup = 3000,
  prior = prior_ma1,
  control = list(adapt_delta=0.95, max_treedepth=15)
)

# save this as rds
saveRDS(fit_ma1, here("Rdata", "fit_ma1.rds"))

fit_ma1 <- readRDS(here("Rdata", "fit_ma1.rds"))

summary(fit_ma1)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: dARR ~ 1 + (1 | p | study_ID) + (1 | gr(es_ID, cov = vcv)) 
##          sigma ~ 1 + (1 | p | study_ID)
##    Data: dat (Number of observations: 1089) 
##   Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup draws = 6000
## 
## Multilevel Hyperparameters:
## ~es_ID (Number of levels: 1089) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## ~study_ID (Number of levels: 150) 
##                                Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                      0.14      0.01     0.11     0.16 1.00
## sd(sigma_Intercept)                0.82      0.07     0.68     0.97 1.00
## cor(Intercept,sigma_Intercept)     0.38      0.12     0.12     0.59 1.00
##                                Bulk_ESS Tail_ESS
## sd(Intercept)                      1103     2257
## sd(sigma_Intercept)                1031     1766
## cor(Intercept,sigma_Intercept)      699     1347
## 
## Regression Coefficients:
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           0.20      0.01     0.17     0.23 1.00     1078     1794
## sigma_Intercept    -2.08      0.09    -2.25    -1.91 1.00     1029     1954
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Note that there is non-zero variance (SD) in the between-study effect on the scale part as well as the location part. Also, there is a positive and significant correaltion (\(\rho_u = 0.38\)) between the location and scale between-study random effects.

3.3.1 Visualzing meta-analytic results

# getting plots
plots_fit_ma1 <- list(
  visualize_fixed_effects(fit_ma1),
  visualize_random_effects(fit_ma1),
  visualize_correlations(fit_ma1)
)


plots_fit_ma1[[1]] / plots_fit_ma1[[2]] / plots_fit_ma1[[3]]

TODO - probably try to get - CV and I2????

3.4 Location-scale meta-regression with a categorical moderator (biological)

#TODO - may need to rerun this model

# biological - meta-regression
form_mr1 <- bf(dARR
            ~ 1  + habitat +
              (1|study_ID) + # this is u_l
              (1|gr(es_ID, cov = vcv)), # this is m
            sigma ~ 1 + habitat
)


# create prior

prior_mr1 <- default_prior(form1, 
                        data = dat, 
                        data2 = list(vcv = vcv),
                        family = gaussian()
)

# fixing the variance to 1 (meta-analysis)
prior_mr1$prior[5] = "constant(1)"
prior_mr1 

# fitting model
fit_mr1 <- brm(form_mr1, 
            data = dat, 
            data2 = list(vcv = vcv),
            chains = 2, 
            cores = 2, 
            iter = 3000, 
            warmup = 2000,
            prior = prior_mr1,
            control = list(adapt_delta = 0.95, max_treedepth = 15)
)

# save this as rds

saveRDS(fit_mr1, here("Rdata", "fit_mr1.rds"))

fit_mr1 <- readRDS(here("Rdata", "fit_mr1.rds"))

summary(fit_mr1)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: dARR ~ 1 + habitat + (1 | study_ID) + (1 | gr(es_ID, cov = vcv)) 
##          sigma ~ 1 + habitat
##    Data: dat (Number of observations: 1089) 
##   Draws: 2 chains, each with iter = 3000; warmup = 2000; thin = 1;
##          total post-warmup draws = 2000
## 
## Multilevel Hyperparameters:
## ~es_ID (Number of levels: 1089) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## ~study_ID (Number of levels: 150) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.13      0.01     0.11     0.16 1.00      543     1180
## 
## Regression Coefficients:
##                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                    0.22      0.02     0.19     0.25 1.00      774
## sigma_Intercept             -1.39      0.03    -1.44    -1.33 1.00     1376
## habitatterrestrial          -0.16      0.04    -0.23    -0.09 1.01      289
## sigma_habitatterrestrial    -1.18      0.08    -1.33    -1.02 1.00     1164
##                          Tail_ESS
## Intercept                    1290
## sigma_Intercept              1678
## habitatterrestrial            781
## sigma_habitatterrestrial     1582
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

3.4.1 Visualzing meta-regression results (biological)

# getting plots
plots_fit_mr1 <- list(
  visualize_fixed_effects(fit_mr1),
  visualize_random_effects(fit_mr1)
)


plots_fit_mr1[[1]] / plots_fit_mr1[[2]]

3.4.2 Visualzing effect sizes with a orchard plot (biological)

# using metafor and use heteroscad model
# see - Nakagawa S, Lagisz M, O'Dea RE, Pottier P, Rutkowska J, Senior AM, Yang Y, Noble DW. orchaRd 2.0: An R package for visualising meta-analyses with orchard plots.
mr1_mod <- rma.mv(yi = dARR, 
                  V =  vcv,
                  mod = ~ habitat,
                  random = list(~1 | study_ID,
                                ~habitat | es_ID),
                  struct = "DIAG",
                  data = dat, 
                  test = "t",
                  sparse = TRUE,
                  control=list(optimizer="optim", optmethod="Nelder-Mead")
)

orchard_plot(mr1_mod, mod = "habitat", group = "study_ID", xlab = "Effect size")

3.5 Location-scale meta-regression with a categorical moderator (methodological)

#TODO - may need to rerun this model

form_mr2 <- bf(dARR
            ~ 1  + method +
              (1|study_ID) + # this is u
              (1|gr(es_ID, cov = vcv)), # this is m
            sigma ~ 1 + method
)


# create prior

prior_mr2 <- default_prior(form_mr2, 
                        data = dat, 
                        data2 = list(vcv = vcv),
                        family = gaussian()
)

# fixing the varaince to 1 (meta-analysis)
prior_mr2$prior[5] = "constant(1)"
prior_mr2 
# fit model

fit_mr2 <- brm(form_mr2, 
            data = dat, 
            data2 = list(vcv = vcv),
            chains = 2, 
            cores = 2, 
            iter = 3000, 
            warmup = 2000,
            #backend = "cmdstanr",
            prior = prior_mr2,
            #threads = threading(9),
            control = list(adapt_delta = 0.99, max_treedepth = 15)
)

# save this as rds

saveRDS(fit1b, here("Rdata", "fit1b.rds"))

fit_mr2 <- readRDS(here("Rdata", "fit_mr2.rds"))

summary(fit_mr2)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: dARR ~ 1 + method + (1 | study_ID) + (1 | gr(es_ID, cov = vcv)) 
##          sigma ~ 1 + method
##    Data: dat (Number of observations: 1089) 
##   Draws: 2 chains, each with iter = 3000; warmup = 2000; thin = 1;
##          total post-warmup draws = 2000
## 
## Multilevel Hyperparameters:
## ~es_ID (Number of levels: 1089) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## ~study_ID (Number of levels: 150) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.12      0.01     0.10     0.15 1.00      718     1373
## 
## Regression Coefficients:
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                   0.24      0.02     0.21     0.28 1.00     1563
## sigma_Intercept            -1.57      0.06    -1.68    -1.45 1.00     1806
## methodpersisitent          -0.07      0.02    -0.10    -0.03 1.00     3050
## sigma_methodpersisitent     0.21      0.07     0.07     0.34 1.00     1967
##                         Tail_ESS
## Intercept                   1667
## sigma_Intercept             1343
## methodpersisitent           1665
## sigma_methodpersisitent     1364
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

3.5.1 Visualzing meta-regression results (methodological)

# getting plots
plots_fit_mr2 <- list(
  visualize_fixed_effects(fit_mr2),
  visualize_random_effects(fit_mr2)
)


plots_fit_mr2[[1]] / plots_fit_mr2[[2]]

3.5.2 Visualzing effect sizes with a orchard plot (methodolgical)

# using metafor and use heteroscad model
# see - Nakagawa S, Lagisz M, O'Dea RE, Pottier P, Rutkowska J, Senior AM, Yang Y, Noble DW. orchaRd 2.0: An R package for visualising meta-analyses with orchard plots.
mr2_mod <- rma.mv(yi = dARR, 
                  V =  vcv,
                  mod = ~ method,
                  random = list(~1 | study_ID,
                                ~method | es_ID),
                  struct = "DIAG",
                  data = dat, 
                  test = "t",
                  sparse = TRUE,
                  control=list(optimizer="optim", optmethod="Nelder-Mead")
)

orchard_plot(mr2_mod, mod = "method", group = "study_ID", xlab = "Effect size")

3.6 Location-scale phylogenetic meta-regression (biological)

TODO - the tree can be obtained from the code ….

3.6.1 Visualzing phylogenetic meta-regression results

4 Example 2: Continuous Moderator (Elevation Dataset)

Data source (illustrative): Midolo, G. et al. (2019). Global patterns of intraspecific leaf trait responses to elevation. Global Change Biology, 25(7), 2485–2498.

4.1 Loading the data

This dataset (elevation.csv) has an effect size (standarised mean difference: SMD or d) and sampling variance Var_dARR. A key moderator is elevation_log (i.e., which elevation on the ln scale, organisms live).

dat <- read.csv(here("data", "elevation.csv"))

dat <- escalc(measure = "ROM", 
              m1i = treatment, 
              m2i = control, 
              sd1i = sd_treatment, 
              sd2i = sd_control, 
              n1i = n_treatment, 
              n2i = n_control, 
              data = dat)


# renaming
dat$study_ID <- as.factor(dat$Study_ID)
dat$es_ID <- as.factor(1:nrow(dat))

# selecting varaibles we need
dat <- dat %>% select(yi, vi, es_ID, study_ID, elevation_log)

#head(dat)
str(dat)
## Classes 'escalc' and 'data.frame':   1294 obs. of  5 variables:
##  $ yi           : num  -0.2071 -0.1881 0.0903 0.0768 -0.0899 ...
##   ..- attr(*, "ni")= int [1:1294] 12 12 12 12 20 20 20 20 90 90 ...
##   ..- attr(*, "measure")= chr "ROM"
##  $ vi           : num  0.022603 0.013296 0.014744 0.016572 0.000585 ...
##  $ es_ID        : Factor w/ 1294 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ study_ID     : Factor w/ 71 levels "1","2","3","4",..: 3 3 3 3 32 32 32 32 36 36 ...
##  $ elevation_log: num  5.52 6.06 4.79 6.06 5.14 ...
##  - attr(*, "yi.names")= chr "yi"
##  - attr(*, "vi.names")= chr "vi"
##  - attr(*, "digits")= Named num [1:9] 4 4 4 4 4 4 4 4 4
##   ..- attr(*, "names")= chr [1:9] "est" "se" "test" "pval" ...

4.2 Creating a variance-covariance matrix for sampling errors

4.2.1 Assuming indepedence among effect sizes

If each effect size has an independent sampling error, we can build a diagonal matrix:

vcv <- diag(dat$vi)
rownames(vcv) <- dat$es_ID
colnames(vcv) <- dat$es_ID

4.3 Multilevel location-scale meta-analysis

#TODO - may need to rerun this model

form_ma2 <- bf(yi
             ~ 1   +
               (1|p|study_ID) + # this is u
               (1|gr(es_ID, cov = vcv)), # this is m
             sigma ~ 1 + (1|p|study_ID)
)



prior_ma2 <- default_prior(form_ma2, 
                         data = dat, 
                         data2 = list(vcv = vcv),
                         family = gaussian()
)

# fixing the varaince to 1 (meta-analysis)
prior_ma2$prior[5] = "constant(1)"
prior_ma2 
# fit model

fit_ma2 <- brm(form_ma2, 
             data = dat, 
             data2 = list(vcv = vcv),
             chains = 2, 
             cores = 2, 
             iter = 6000, 
             warmup = 3000,
             prior = prior_ma2,
             control = list(adapt_delta = 0.95, max_treedepth = 15)
)

# save this as rds

saveRDS(fit_ma2, here("Rdata", "fit_ma2.rds"))

fit_ma2 <- readRDS(here("Rdata", "fit_ma2.rds"))

summary(fit_ma2)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: yi ~ 1 + (1 | p | study_ID) + (1 | gr(es_ID, cov = vcv)) 
##          sigma ~ 1 + (1 | p | study_ID)
##    Data: dat (Number of observations: 1294) 
##   Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup draws = 6000
## 
## Multilevel Hyperparameters:
## ~es_ID (Number of levels: 1294) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## ~study_ID (Number of levels: 71) 
##                                Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                      0.11      0.02     0.08     0.14 1.00
## sd(sigma_Intercept)                0.74      0.09     0.60     0.93 1.00
## cor(Intercept,sigma_Intercept)     0.39      0.15     0.07     0.65 1.00
##                                Bulk_ESS Tail_ESS
## sd(Intercept)                      1170     1874
## sd(sigma_Intercept)                1091     1415
## cor(Intercept,sigma_Intercept)      667     1379
## 
## Regression Coefficients:
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           0.02      0.02    -0.01     0.06 1.00      746     1433
## sigma_Intercept    -1.84      0.10    -2.05    -1.64 1.01      960     1844
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

4.3.1 Visualzing meta-analytic results

# getting plots
plots_fit_ma2 <- list(
  visualize_fixed_effects(fit_ma2),
  visualize_random_effects(fit_ma2),
  visualize_correlations(fit_ma2)
)


plots_fit_ma2[[1]] / plots_fit_ma2[[2]] / plots_fit_ma2[[3]]

TODO - probably try to get - CV and I2????

4.4 Location-scale meta-regression with a contnious moderator

#TODO - may need to rerun this model

# SMD = d 
form_mr3 <- bf(yi
            ~ 1  + elevation_log +
              (1|study_ID) + # this is u
              (1|gr(es_ID, cov = vcv)), # this is m
            sigma ~ 1 + elevation_log
)


# create prior


prior_mr3 <- default_prior(form_mr3, 
                        data = dat, 
                        data2 = list(vcv = vcv),
                        family = gaussian()
)

# fixing the varaince to 1 (meta-analysis)
prior_mr3$prior[5] = "constant(1)"
prior_mr3 

# fit model
fit_mr3 <- brm(form_mr3, 
            data = dat, 
            data2 = list(vcv = vcv),
            chains = 2, 
            cores = 2, 
            iter = 6000, 
            warmup = 3000,
            prior = prior_mr3,
            control = list(adapt_delta = 0.95, max_treedepth = 15)
)

# save this as rds

saveRDS(fit_mr3, here("Rdata", "fit_mr3.rds"))

fit_mr3 <- readRDS(here("Rdata", "fit_mr3.rds"))

summary(fit_mr3)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: yi ~ 1 + elevation_log + (1 | study_ID) + (1 | gr(es_ID, cov = vcv)) 
##          sigma ~ 1 + elevation_log
##    Data: dat (Number of observations: 1294) 
##   Draws: 2 chains, each with iter = 6000; warmup = 3000; thin = 1;
##          total post-warmup draws = 6000
## 
## Multilevel Hyperparameters:
## ~es_ID (Number of levels: 1294) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## ~study_ID (Number of levels: 71) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.13      0.02     0.10     0.16 1.00     1104     2260
## 
## Regression Coefficients:
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept              -0.13      0.05    -0.23    -0.02 1.00     2213     3632
## sigma_Intercept        -3.73      0.18    -4.09    -3.39 1.00     3566     4319
## elevation_log           0.03      0.01     0.01     0.04 1.00     3332     4231
## sigma_elevation_log     0.36      0.03     0.30     0.42 1.00     3744     4498
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

4.4.1 Visualzing meta-regression results

# getting plots
plots_fit_mr3 <- list(
  visualize_fixed_effects(fit_mr3),
  visualize_random_effects(fit_mr3)
)


plots_fit_mr3[[1]] / plots_fit_mr3[[2]]

4.4.2 Visualzing effect sizes with a bubble plot

# using metafor and orchaRd
# see - Nakagawa S, Lagisz M, O'Dea RE, Pottier P, Rutkowska J, Senior AM, Yang Y, Noble DW. orchaRd 2.0: An R package for visualising meta-analyses with orchard plots.
mr3_mod <- rma.mv(yi = yi, 
                  V =  vcv,
                  mod = ~ elevation_log,
                  random = list(~1 | study_ID,
                                ~1 | es_ID),
                  #struct = "DIAG",
                  data = dat, 
                  test = "t",
                  sparse = TRUE,
                  control=list(optimizer="optim", optmethod="Nelder-Mead")
)

bubble_plot(mr3_mod, mod = "elevation_log", group = "study_ID", g = TRUE, 
            xlab = "ln(elevation) (moderator)",
            ylab = "SDM (effect size)")

One can clearly see an increase in variation as evvelation_log goes up.

5 Example 3: Publication Bias (Small-Study, Decline, Proteus)

Data source (illustrative): Neuschulz, E. L. et al. (2016). Pollination and seed dispersal are the most threatened processes of plant regeneration. Scientific Reports, 6, 29839.

We illustrate how to incorporate se (or ()) and cyear (centered publication year) in both the location and scale parts:

6 Loading the data

This dataset (seed.csv) has an effect size (SND) and sampling variance. Two key moderators are: (1) se (standard error for sampling error to model the small-study effect and divergence) and cyear (centred publication year to model the decline effect and Proteus effect).

dat <- read.csv(here("data", "seed.csv"))

dat$study_ID <- as.factor(dat$study)
dat$es_ID <- as.factor(1:nrow(dat))
dat$se <- sqrt(dat$var.eff.size) # SE (sampling error SD)
dat$smd <- dat$eff.size # effect isze
dat$cyear <- scale(dat$study.year, center = TRUE, scale = FALSE) # centred year


# selecting varaibles we need
dat <- dat %>% select(smd, var.eff.size, es_ID, study_ID, se, cyear)

#head(dat)
str(dat)
## 'data.frame':    98 obs. of  6 variables:
##  $ smd         : num  -0.65 2.24 -1.39 -1.08 0.42 0.34 -0.05 -5.53 -1.07 -0.46 ...
##  $ var.eff.size: num  0.6 0.23 0.17 0.16 0.14 0.14 0.13 0.52 0.49 0.42 ...
##  $ es_ID       : Factor w/ 98 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ study_ID    : Factor w/ 40 levels "Albrecht et al. 2012",..: 34 20 20 20 20 20 20 33 40 40 ...
##  $ se          : num  0.775 0.48 0.412 0.4 0.374 ...
##  $ cyear       : num [1:98, 1] -13.59 -7.59 -7.59 -7.59 -7.59 ...
##   ..- attr(*, "scaled:center")= num 2008

6.1 Creating a variance-covariance matrix for sampling errors

6.1.1 Assuming indepedence among effect sizes

If each effect size has an independent sampling error, we can build a diagonal matrix:

vcv <- diag(dat$var.eff.size)
rownames(vcv) <- dat$es_ID
colnames(vcv) <- dat$es_ID

6.2 Multilevel location-scale meta-analysis with correlation

#TODO - may need to rerun this model

form_ma3 <- bf(smd
            ~ 1   +
              (1|study_ID) + # this is u_l
              (1|gr(es_ID, cov = vcv)), # this is m
            sigma ~ 1 + (1|study_ID)
)



prior_ma3 <- default_prior(form_ma3, 
                        data = dat, 
                        data2 = list(vcv = vcv),
                        family = gaussian()
)


# fixing the varaince to 1 (meta-analysis)
prior_ma3$prior[3] = "constant(1)"
prior_ma3

# fit model
fit_ma3 <- brm(form_ma3, 
            data = dat, 
            data2 = list(vcv = vcv),
            chains = 2, 
            cores = 2, 
            iter = 30000, 
            warmup = 5000,
            prior = prior_ma3,
            control = list(adapt_delta = 0.99, max_treedepth = 15)
)

summary(fit_ma3)

# save this as rds

saveRDS(fit_ma3, here("Rdata", "fit_ma3.rds"))

fit_ma3 <- readRDS(here("Rdata", "fit_ma3.rds"))

summary(fit_ma3)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: smd ~ 1 + (1 | p | study_ID) + (1 | gr(es_ID, cov = vcv)) 
##          sigma ~ 1 + (1 | p | study_ID)
##    Data: dat (Number of observations: 98) 
##   Draws: 2 chains, each with iter = 35000; warmup = 5000; thin = 1;
##          total post-warmup draws = 60000
## 
## Multilevel Hyperparameters:
## ~es_ID (Number of levels: 98) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.63      0.23     1.10     1.85 2.09        4       54
## 
## ~study_ID (Number of levels: 40) 
##                                Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                      0.58      0.15     0.36     0.96 1.69
## sd(sigma_Intercept)                1.97      0.82     0.18     2.78 1.91
## cor(Intercept,sigma_Intercept)     0.46      0.61    -0.87     0.98 2.01
##                                Bulk_ESS Tail_ESS
## sd(Intercept)                         8       39
## sd(sigma_Intercept)                   3       25
## cor(Intercept,sigma_Intercept)        3       15
## 
## Regression Coefficients:
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          -0.35      0.13    -0.70    -0.15 1.74       15       37
## sigma_Intercept    -4.52      2.21    -6.78    -0.92 2.03        3       36
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

The parameters are not mixing well and we have covergence issues which are also clear from the fig below.

6.2.1 Visualzing meta-analytic results (without correlation)

# getting plots
plots_fit_ma3 <- list(
  visualize_fixed_effects(fit_ma3),
  visualize_random_effects(fit_ma3),
  visualize_correlations(fit_ma3)
)


plots_fit_ma3[[1]] / plots_fit_ma3[[2]] / plots_fit_ma3[[3]]

6.3 Multilevel location-scale meta-analysis without the between-study corrleation

#TODO - may need to rerun this model

# no correlation 
# this runs but the one with correlation does not mix well

form_ma4 <- bf(smd
            ~ 1   +
              (1|study_ID) + # this is u
              (1|gr(es_ID, cov = vcv)), # this is m
            sigma ~ 1 + (1|study_ID)
)



prior_ma4 <- default_prior(form_ma4, 
                        data = dat, 
                        data2 = list(vcv = vcv),
                        family = gaussian()
)


# fixing the varaince to 1 (meta-analysis)
prior_ma4$prior[3] = "constant(1)"
prior_ma4 
# fit model

fit_ma4 <- brm(form_ma4, 
            data = dat, 
            data2 = list(vcv = vcv),
            chains = 2, 
            cores = 2, 
            iter = 30000, 
            warmup = 5000,
            #backend = "cmdstanr",
            prior = prior_ma4,
            #threads = threading(9),
            control = list(adapt_delta = 0.99, max_treedepth = 15)
)

summary(fit_ma4)

# save this as rds

saveRDS(fit4, here("Rdata", "fit_ma4.rds"))

fit_ma4 <- readRDS(here("Rdata", "fit_ma4.rds"))

summary(fit_ma4)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: smd ~ 1 + (1 | study_ID) + (1 | gr(es_ID, cov = vcv)) 
##          sigma ~ 1 + (1 | study_ID)
##    Data: dat (Number of observations: 98) 
##   Draws: 2 chains, each with iter = 30000; warmup = 5000; thin = 1;
##          total post-warmup draws = 50000
## 
## Multilevel Hyperparameters:
## ~es_ID (Number of levels: 98) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## ~study_ID (Number of levels: 40) 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)           0.66      0.15     0.41     0.98 1.00     2164      816
## sd(sigma_Intercept)     1.53      0.51     0.72     2.72 1.02      387      277
## 
## Regression Coefficients:
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept          -0.46      0.15    -0.77    -0.17 1.00     2895     5918
## sigma_Intercept    -1.40      0.66    -2.97    -0.40 1.02      374      189
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

The parameters are not mixing well and we have covergence issues which are also clear from the fig below.

6.3.1 Visualzing meta-analytic results (without correlation)

# getting plots
plots_fit_ma4 <- list(
  visualize_fixed_effects(fit_ma4),
  visualize_random_effects(fit_ma4),
  visualize_correlations(fit_ma4)
)


plots_fit_ma4[[1]] / plots_fit_ma4[[2]] 

TODO - it is mixing well

6.4 Location-scale meta-regression to model publication bias

#TODO - may need to rerun this model

form_mr4 <- bf(smd
             ~ 1   +
               (1|p|study_ID) + # this is u
               (1|gr(es_ID, cov = vcv)), # this is m
             sigma ~ 1 + (1|p|study_ID)
)



prior_mr4 <- default_prior(form_mr4, 
                         data = dat, 
                         data2 = list(vcv = vcv),
                         family = gaussian()
)

# fixing the varaince to 1 (meta-analysis)
prior4b$prior[5] = "constant(1)"
prior4b 
# fit model

fit_mr4 <- brm(form_mr4, 
             data = dat, 
             data2 = list(vcv = vcv),
             chains = 2, 
             cores = 2, 
             iter =35000, 
             warmup = 5000,
             #backend = "cmdstanr",
             prior = prior_mr4,
             #threads = threading(9),
             control = list(adapt_delta = 0.99, max_treedepth = 15)
)

summary(fit_mr4)

# save this as rds

saveRDS(fit_mr4, here("Rdata", "fit_mr4.rds"))

fit_mr4 <- readRDS(here("Rdata", "fit_mr4.rds"))

summary(fit_mr4)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: smd ~ 1 + cyear + se + (1 | study_ID) + (1 | gr(es_ID, cov = vcv)) 
##          sigma ~ 1 + cyear + se
##    Data: dat (Number of observations: 98) 
##   Draws: 2 chains, each with iter = 35000; warmup = 5000; thin = 1;
##          total post-warmup draws = 60000
## 
## Multilevel Hyperparameters:
## ~es_ID (Number of levels: 98) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## ~study_ID (Number of levels: 40) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.66      0.13     0.44     0.95 1.00    15622    27363
## 
## Regression Coefficients:
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept           0.10      0.33    -0.53     0.75 1.00    19025    27181
## sigma_Intercept    -2.30      0.61    -3.59    -1.19 1.00     1933     2693
## cyear              -0.04      0.04    -0.12     0.04 1.00    28664    36930
## se                 -0.89      0.58    -2.06     0.23 1.00    15695    22795
## sigma_cyear        -0.19      0.06    -0.31    -0.09 1.00     1317     1709
## sigma_se            2.13      0.67     0.74     3.41 1.00     4982     6298
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

6.4.1 Visualzing meta-regression results

# getting plots
plots_fit_mr4 <- list(
  visualize_fixed_effects(fit_mr4),
  visualize_random_effects(fit_mr4)
)


plots_fit_mr4[[1]] / plots_fit_mr4[[2]]

6.4.2 Visualzing effect sizes with a bubble plot

# using metafor and orchaRd
# see - Nakagawa S, Lagisz M, O'Dea RE, Pottier P, Rutkowska J, Senior AM, Yang Y, Noble DW. orchaRd 2.0: An R package for visualising meta-analyses with orchard plots.
mr4_mod <- rma.mv(yi = smd, 
                  V =  vcv,
                  mod = ~ se + cyear,
                  random = list(~1 | study_ID,
                                ~1 | es_ID),
                  #struct = "DIAG",
                  data = dat, 
                  test = "t",
                  sparse = TRUE,
                  control=list(optimizer="optim", optmethod="Nelder-Mead")
)

bubble_plot(mr4_mod, mod = "se", group = "study_ID", g = TRUE, 
            xlab = "Standard error (SE: moderator)",
            ylab = "SDM (effect size)")

bubble_plot(mr4_mod, mod = "cyear", group = "study_ID", g = TRUE, 
            xlab = "centered publicaiton year (moderator)",
            ylab = "SDM (effect size)")

One …..

7 Comparing with other packages (metafor and blsmeta)

Here we compare results of three packages. Note the followings:

• metafor can only ran the random-effect meta-regression
• blsmeta can run multilevel but only to the 3 levels but it cannot take random effects on the scale part. Instead, blsmeta can have regression formulas for each of variance components (apart from the sampling error variance; this is a different model and conceptualization from our proposed model). For the details reading the following two papers

TODO - just cut and pasted my test2.R file - need to think of some subheadings

#####################
# random effect model
#####################

# metafor 1
#------------
ma1 <- rma(yi = dARR, 
           vi = Var_dARR, 
           test="t",
           data = dat)

summary(ma1)


# metafor 2
#------------
ma2 <- rma.mv(yi = dARR, 
              V = vcv,
              random = ~ 1 | es_ID,
              test="t",
              data = dat)

summary(ma2)


# brms 1
#------------

form3 <- bf(dARR | se(si) ~ 1 +  (1|es_ID) # this is e
)

prior3 <- default_prior(form3, 
                        data = dat, 
                        #data2 = list(vcv = vcv),
                        family = gaussian()
)

ma3 <- brm(form3, 
            data = dat, 
            #data2 = list(vcv = vcv),
            chains = 2, 
            cores = 2, 
            iter = 35000, 
            warmup = 5000,
            #backend = "cmdstanr",
            prior = prior3,
            #threads = threading(9),
            control = list(adapt_delta = 0.99, max_treedepth = 20)
)


# save as rds
saveRDS(ma3, here("Rdata", "ma3.rds"))

# load rds
ma3 <- readRDS(here("Rdata", "ma3.rds"))


print(summary(ma3), digits = 4)

# brms 2
#------------
form4 <- bf(dARR  ~ 1 +  (1|es_ID) # this is e
            + fcor(vcv)
)

prior4 <- default_prior(form4, 
                        data = dat, 
                        data2 = list(vcv = vcv),
                        family = gaussian()
)


# we need to fix the variance to 1 = fcor(vcv)
prior4$prior[5] = "constant(1)"
prior4 

ma4 <- brm(form4, 
           data = dat, 
           data2 = list(vcv = vcv),
           chains = 2, 
           cores = 2, 
           iter = 15000, 
           warmup = 5000,
           #backend = "cmdstanr",
           prior = prior4,
           #threads = threading(9),
           #control = list(adapt_delta = 0.95, max_treedepth = 15)
)

# save as rds
saveRDS(ma4, here("Rdata", "ma4.rds"))

# load rds
ma4 <- readRDS(here("Rdata", "ma4.rds"))

print(summary(ma4), digit = 4)


# brm 3 (using Shinichi's trick)
#------------
form5 <- bf(dARR  ~ 1 +   (1|gr(es_ID, cov = vcv)), # this is m
            # residual will be es_ID SD
)

prior5 <- default_prior(form5, 
                        data = dat, 
                        data2 = list(vcv = vcv),
                        family = gaussian()
)

# fixing the varaince to 1 (meta-analysis)
prior5$prior[3] = "constant(1)"
prior5 

ma5 <- brm(form5, 
           data = dat, 
           data2 = list(vcv = vcv),
           chains = 2, 
           cores = 2, 
           iter = 35000, 
           warmup = 5000,
           prior = prior5
)

print(summary(ma5), digits = 4)


# above shows my trick works - let's move to l-s meta-regressoin

# metafor
#------------
mr1 <- rma(yi = dARR, 
            vi = Var_dARR, 
            mods = ~ habitat,
            scale = ~ habitat,
            test="t",
            data = dat)

summary(mr1)


# brms
#------------
form2 <- bf(dARR
            ~ 1   + habitat + 
              (1|gr(es_ID, cov = vcv)), # this is m
            sigma ~ 1 + habitat
)


prior2 <- default_prior(form2, 
                        data = dat, 
                        data2 = list(vcv = vcv),
                        family = gaussian()
)


# fixing the varaince to 1 (meta-analysis)
prior2$prior[5] = "constant(1)"
prior2 
# fit model

mr2 <- brm(form2, 
            data = dat, 
            data2 = list(vcv = vcv),
            chains = 2, 
            cores = 2, 
            iter = 35000, 
            warmup = 5000,
            #backend = "cmdstanr",
            prior = prior2,
            #threads = threading(9),
            control = list(adapt_delta = 0.95, max_treedepth = 15)
)

# save as rds

saveRDS(mr2, here("Rdata", "mr2.rds"))

# load the rds

mr2 <- readRDS(here("Rdata", "mr2.rds"))

summary(mr2)

# blmata


res0 <- blsmeta(yi = dARR,
               vi = Var_dARR,
               es_id = es_ID,
               study_id = study_ID,
               mods_scale2 = ~ 1,
               data = dat)


###############
# common mistake - the model which does not run

# brms

form1 <- bf(dARR | se(si) ~ 1 + habitat, # this is e
            sigma ~ 1 + habitat
)

prior1 <- default_prior(form1, 
                        data = dat, 
                        #data2 = list(vcv = vcv),
                        family = gaussian()
)

ma1 <- brm(form1, 
           data = dat, 
           #data2 = list(vcv = vcv),
           chains = 2, 
           cores = 2, 
           iter = 35000, 
           warmup = 5000,
           #backend = "cmdstanr",
           prior = prior1,
           #threads = threading(9),
           control = list(adapt_delta = 0.99, max_treedepth = 20)
)

#############
## 3 levels
################

# biological - meta-regression
form1 <- bf(dARR
            ~ 1  + habitat +
              (1|study_ID) + # this is u
              (1|gr(es_ID, cov = vcv)), # this is m
            sigma ~ 1 + habitat
)


# create prior

prior1 <- default_prior(form1, 
                        data = dat, 
                        data2 = list(vcv = vcv),
                        family = gaussian()
)

# fixing the varaince to 1 (meta-analysis)
prior1$prior[5] = "constant(1)"
prior1 
# fit model

res1 <- brm(form1, 
            data = dat, 
            data2 = list(vcv = vcv),
            chains = 2, 
            cores = 2, 
            iter = 3000, 
            warmup = 2000,
            #backend = "cmdstanr",
            prior = prior1,
            #threads = threading(9),
            control = list(adapt_delta = 0.95, max_treedepth = 15)
)

# save this as rds

saveRDS(res1, here("Rdata", "res1.rds"))

# read in rds
res1_brms <- readRDS(here("Rdata", "res1.rds"))


summary(res1_brms)

# blmeta

res1a <- blsmeta(yi = dARR,
               vi = Var_dARR,
               es_id = es_ID,
               study_id = study_ID,
               mods = ~ 1 + habitat,
               mods_scale2 = ~ 1 + habitat,
               data = dat)
#save

saveRDS(res1a, here("Rdata", "res1a.rds"))

# read in rds

res1a <- readRDS(here("Rdata", "res1a.rds"))

print(res1a)


# try method

dat$method <- ifelse(dat$exp_design == c("A", "B", "C"), "initial", "persisitent")



form1b <- bf(dARR
             ~ 1  + method +
               (1|study_ID) + # this is u
               (1|gr(es_ID, cov = vcv)), # this is m
             sigma ~ 1 + method
)


# create prior

prior1b <- default_prior(form1b, 
                         data = dat, 
                         data2 = list(vcv = vcv),
                         family = gaussian()
)

# fixing the varaince to 1 (meta-analysis)
prior1b$prior[5] = "constant(1)"
prior1b 
# fit model

res2 <- brm(form1b, 
             data = dat, 
             data2 = list(vcv = vcv),
             chains = 2, 
             cores = 2, 
             iter = 3000, 
             warmup = 2000,
             #backend = "cmdstanr",
             prior = prior1b,
             #threads = threading(9),
             control = list(adapt_delta = 0.99, max_treedepth = 15)
)

summary(res2)

# save this as rds
saveRDS(res2, here("Rdata", "res2.rds"))
# reading rsd

res2 <- readRDS(here("Rdata", "res2.rds"))

summary(res2)

# blsmeta


res2a <- blsmeta(yi = dARR,
                 vi = Var_dARR,
                 es_id = es_ID,
                 study_id = study_ID,
                 mods = ~ 1 + method,
                 mods_scale2 = ~ 1 + method,
                 data = dat)
#save

#save

saveRDS(res2a, here("Rdata", "res2a.rds"))

# read in rds

res2a <- readRDS(here("Rdata", "res2a.rds"))

print(res2a)

Note

xxxxx

It’s important…

Important

In …

8 Software and package versions

sessionInfo() %>% pander()

R version 4.4.2 (2024-10-31)

Platform: aarch64-apple-darwin20

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: blsmeta(v.0.1.0), metafor(v.4.6-0), numDeriv(v.2016.8-1.1), metadat(v.1.2-0), Matrix(v.1.7-0), brms(v.2.21.0), Rcpp(v.1.0.12), pander(v.0.6.5), orchaRd(v.2.0), patchwork(v.1.2.0), here(v.1.0.1), tidybayes(v.3.0.7), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.2), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.2.1), ggplot2(v.3.5.1), tidyverse(v.2.0.0) and pacman(v.0.5.1)

loaded via a namespace (and not attached): gridExtra(v.2.3), inline(v.0.3.19), sandwich(v.3.1-0), rlang(v.1.1.4), magrittr(v.2.0.3), multcomp(v.1.4-25), matrixStats(v.1.4.1), compiler(v.4.4.2), mgcv(v.1.9-1), loo(v.2.8.0), vctrs(v.0.6.5), reshape2(v.1.4.4), bayeslincom(v.1.3.0), pkgconfig(v.2.0.3), arrayhelpers(v.1.1-0), fastmap(v.1.2.0), backports(v.1.5.0), labeling(v.0.4.3), utf8(v.1.2.4), rmarkdown(v.2.27), tzdb(v.0.4.0), ggbeeswarm(v.0.7.2), xfun(v.0.45), jsonlite(v.1.8.8), parallel(v.4.4.2), R6(v.2.5.1), stringi(v.1.8.4), StanHeaders(v.2.32.10), rjags(v.4-16), estimability(v.1.5.1), rstan(v.2.32.6), knitr(v.1.48), zoo(v.1.8-12), bayesplot(v.1.11.1), splines(v.4.4.2), timechange(v.0.3.0), tidyselect(v.1.2.1), rstudioapi(v.0.16.0), abind(v.1.4-5), yaml(v.2.3.9), codetools(v.0.2-20), pkgbuild(v.1.4.4), lattice(v.0.22-6), plyr(v.1.8.9), withr(v.3.0.0), bridgesampling(v.1.1-2), posterior(v.1.6.0), coda(v.0.19-4.1), evaluate(v.0.24.0), survival(v.3.7-0), RcppParallel(v.5.1.9), ggdist(v.3.3.2), pillar(v.1.9.0), tensorA(v.0.36.2.1), checkmate(v.2.3.2), stats4(v.4.4.2), distributional(v.0.5.0), generics(v.0.1.3), rprojroot(v.2.0.4), mathjaxr(v.1.6-0), hms(v.1.1.3), rstantools(v.2.4.0), munsell(v.0.5.1), scales(v.1.3.0), xtable(v.1.8-4), glue(v.1.7.0), emmeans(v.1.10.4), tools(v.4.4.2), mvtnorm(v.1.2-5), grid(v.4.4.2), QuickJSR(v.1.3.1), colorspace(v.2.1-0), nlme(v.3.1-165), beeswarm(v.0.4.0), vipor(v.0.4.7), latex2exp(v.0.9.6), cli(v.3.6.3), fansi(v.1.0.6), svUnit(v.1.0.6), Brobdingnag(v.1.2-9), gtable(v.0.3.5), digest(v.0.6.36), TH.data(v.1.1-2), htmlwidgets(v.1.6.4), farver(v.2.1.2), htmltools(v.0.5.8.1), lifecycle(v.1.0.4) and MASS(v.7.3-61)